Brian E. J. Rose, University at Albany
This document uses the interactive IPython notebook
format (now also called Jupyter
). The notes can be accessed in several different ways:
github
at https://github.com/brian-rose/ClimateModeling_coursewareMany of these notes make use of the climlab
package, available at https://github.com/brian-rose/climlab
The surface of the Earth is the boundary between the atmosphere and the land, ocean, or ice. Understanding the energy fluxes across the surface are very important for three main reasons:
The energy budget at the surface is more complex that the budget at the top of the atmosphere. At the TOA the only energy transfer mechanisms are radiative (shortwave and longwave). At the surface, in addition to radiation we need to consider fluxes of energy by conduction and by convection of heat and moisture through turbulent fluid motion.
We will denote the net upward energy flux at the surface as $F_S$.
As we mentioned back in Lecture 13 on heat transport, there are four principal contributions to $F_S$:
Wherever $F_S \ne 0$, there is a net flux of energy between the atmosphere and the surface below. This implies either that there is heat storage / release occuring below the surface (e.g. warming or cooling of water, melting of snow and ice), and/or there is horizontal heat transport by fluid motions occuring below the surface (ocean circulation, groundwater flow).
All of these terms are small globally but can be significant locally or seasonally.
We will examine the surface budget in the CESM slab ocean simulations. The advantage of looking at surface fluxes in a model rather than observations is that the model fluxes are completely consistent with the model climate, so that the net flux $F_S$ will be a meaningful measure of the heat storage in the system.
The model also gives us an opportunity to look at how the surface budget reponds to global warming under a doubling of CO$_2$.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
from climlab import constants as const
datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/"
endstr = "/entry.das"
topo = nc.Dataset( datapath + 'som_input/USGS-gtopo30_1.9x2.5_remap_c050602.nc' + endstr )
runlist = ['control', '2xCO2']
runs = {}
for run in runlist:
runstr = 'som_' + run
path = datapath + runstr + '/' + runstr + '.cam.h0.clim.nc' + endstr
runs[run] = nc.Dataset( path )
lat = runs['control'].variables['lat'][:]
lon = runs['control'].variables['lon'][:]
lev = runs['control'].variables['lev'][:]
# annual mean surface energy budget terms, all defined as positive up (from ocean to atmosphere)
surface_budget = {}
for name, run in runs.iteritems():
budget = {}
budget['LHF'] = np.mean(run.variables['LHFLX'][:], axis=0)
budget['SHF'] = np.mean(run.variables['SHFLX'][:], axis=0)
budget['LWsfc'] = np.mean(run.variables['FLNS'][:], axis=0)
budget['LWsfc_clr'] = np.mean(run.variables['FLNSC'][:], axis=0)
budget['SWsfc'] = -np.mean(run.variables['FSNS'][:], axis=0)
budget['SWsfc_clr'] = -np.mean(run.variables['FSNSC'][:], axis=0)
budget['SnowFlux'] = (np.mean(run.variables['PRECSC'][:]+run.variables['PRECSL'][:], axis=0)
*const.rho_w*const.Lhfus)
budget['NetRad'] = budget['LWsfc'] + budget['SWsfc'] # net upward radiation from surface
# net upward surface heat flux
budget['Net'] = (budget['NetRad'] + budget['LHF'] +
budget['SHF'] + budget['SnowFlux'])
budget['LWsfc_cld'] = budget['LWsfc'] - budget['LWsfc_clr'] # all the surface radiation terms are defined positive up
budget['SWsfc_cld'] = budget['SWsfc'] - budget['SWsfc_clr']
budget['NetRad_clr'] = budget['LWsfc_clr'] + budget['SWsfc_clr']
budget['NetRad_cld'] = budget['NetRad'] - budget['NetRad_clr']
surface_budget[name] = budget
# Compute anomalies for all terms
anom = {}
for name in surface_budget['control']:
anom[name] = surface_budget['2xCO2'][name] - surface_budget['control'][name]
surface_budget['anom'] = anom
# Also compute zonal averages
zonal_budget = {}
for run, budget in surface_budget.iteritems():
thiszon = {}
for field in budget:
thiszon[field] = np.mean(budget[field], axis=1)
zonal_budget[run] = thiszon
fig = plt.figure(figsize=(16,5))
ax1 = fig.add_subplot(1, 2, 1)
cax1 = ax1.pcolormesh(lon, lat, surface_budget['control']['Net'],
cmap=plt.cm.seismic, vmin=-200., vmax=200. )
fig.colorbar(cax1)
ax1.set_title('Annual mean net surface heat flux (+ up) - CESM control')
ax2 = fig.add_subplot(1, 2, 2)
cax2 = ax2.pcolormesh(lon, lat, surface_budget['anom']['Net'],
cmap=plt.cm.seismic, vmin=-20., vmax=20. )
fig.colorbar(cax2)
ax2.set_title('Anomaly after CO2 doubling')
for ax in [ax1, ax2]:
ax.set_xlim(0, 360)
ax.set_ylim(-90, 90)
ax.contour( lon, lat, topo.variables['LANDFRAC'][:], [0.5,0.5], colors='k');
These figures show the annual mean net upward flux $F_S$ in the control run, and the anomaly after greenhouse warming.
Some notable points about the control state:
After greenhouse warming:
fieldlist = ['SWsfc', 'LWsfc', 'LHF', 'SHF', 'Net']
fig = plt.figure(figsize=(16,5))
ax1 = fig.add_subplot(1, 2, 1)
for field in fieldlist:
ax1.plot(lat, zonal_budget['control'][field], label=field)
ax1.set_title('Components of ANNUAL surface energy budget (+ up) - CESM control')
ax2 = fig.add_subplot(1, 2, 2)
for field in fieldlist:
ax2.plot(lat, zonal_budget['anom'][field], label=field)
ax2.set_title('Anomaly after CO2 doubling')
for ax in [ax1, ax2]:
ax.set_xlim(-90, 90)
ax.grid()
ax.legend(loc='lower right')
In these graphs, the magenta curve labeled "Net" is the net flux $F_S$. It is just the zonal average of the maps from the previous figure, and shows the ocean heat uptake at the equator and release in mid- to high latitudes.
More interestingly, these graphs show the contribution of the various terms to $F_S$. They are all plotted as positive up. A negative value thus indicates heating of the surface, and a positive value indicates a cooling of the surface.
Key points about the control simulation:
After greenhouse warming
We will compute the budgets for the months of January and July, and plot their differences.
jan_budget = {}
for name, run in runs.iteritems():
budget = {}
budget['LHF'] = run.variables['LHFLX'][0,:,:]
budget['SHF'] = run.variables['SHFLX'][0,:,:]
budget['LWsfc'] = run.variables['FLNS'][0,:,:]
budget['LWsfc_clr'] = run.variables['FLNSC'][0,:,:]
budget['SWsfc'] = -run.variables['FSNS'][0,:,:]
budget['SWsfc_clr'] = -run.variables['FSNSC'][0,:,:]
budget['SnowFlux'] = ((run.variables['PRECSC'][0,:,:]+run.variables['PRECSL'][0,:,:])
*const.rho_w*const.Lhfus)
budget['NetRad'] = budget['LWsfc'] + budget['SWsfc'] # net upward radiation from surface
# net upward surface heat flux
budget['Net'] = (budget['NetRad'] + budget['LHF'] +
budget['SHF'] + budget['SnowFlux'])
budget['LWsfc_cld'] = budget['LWsfc'] - budget['LWsfc_clr'] # all the surface radiation terms are defined positive up
budget['SWsfc_cld'] = budget['SWsfc'] - budget['SWsfc_clr']
budget['NetRad_clr'] = budget['LWsfc_clr'] + budget['SWsfc_clr']
budget['NetRad_cld'] = budget['NetRad'] - budget['NetRad_clr']
jan_budget[name] = budget
jul_budget = {}
for name, run in runs.iteritems():
budget = {}
budget['LHF'] = run.variables['LHFLX'][6,:,:]
budget['SHF'] = run.variables['SHFLX'][6,:,:]
budget['LWsfc'] = run.variables['FLNS'][6,:,:]
budget['LWsfc_clr'] = run.variables['FLNSC'][6,:,:]
budget['SWsfc'] = -run.variables['FSNS'][6,:,:]
budget['SWsfc_clr'] = -run.variables['FSNSC'][6,:,:]
budget['SnowFlux'] = ((run.variables['PRECSC'][6,:,:]+run.variables['PRECSL'][6,:,:])
*const.rho_w*const.Lhfus)
budget['NetRad'] = budget['LWsfc'] + budget['SWsfc'] # net upward radiation from surface
# net upward surface heat flux
budget['Net'] = (budget['NetRad'] + budget['LHF'] +
budget['SHF'] + budget['SnowFlux'])
budget['LWsfc_cld'] = budget['LWsfc'] - budget['LWsfc_clr'] # all the surface radiation terms are defined positive up
budget['SWsfc_cld'] = budget['SWsfc'] - budget['SWsfc_clr']
budget['NetRad_clr'] = budget['LWsfc_clr'] + budget['SWsfc_clr']
budget['NetRad_cld'] = budget['NetRad'] - budget['NetRad_clr']
jul_budget[name] = budget
# July minus January
julminusjan_budget = {}
for run in runs:
anom = {}
for name in jan_budget[run]:
anom[name] = jul_budget[run][name] - jan_budget[run][name]
julminusjan_budget[run] = anom
# Compute anomalies for all terms
anom = {}
for name in julminusjan_budget['control']:
anom[name] = julminusjan_budget['2xCO2'][name] - julminusjan_budget['control'][name]
julminusjan_budget['anom'] = anom
# Also compute zonal averages
julminusjan_zonal_budget = {}
for run, budget in julminusjan_budget.iteritems():
thiszon = {}
for field in budget:
thiszon[field] = np.mean(budget[field], axis=1)
julminusjan_zonal_budget[run] = thiszon
fieldlist = ['SWsfc', 'LWsfc', 'LHF', 'SHF', 'Net']
fig = plt.figure(figsize=(16,5))
ax1 = fig.add_subplot(1, 2, 1)
for field in fieldlist:
ax1.plot(lat, julminusjan_zonal_budget['control'][field], label=field)
ax1.set_title('Components of JUL-JAN surface energy budget (+ up) - CESM control')
ax2 = fig.add_subplot(1, 2, 2)
for field in fieldlist:
ax2.plot(lat, julminusjan_zonal_budget['anom'][field], label=field)
ax2.set_title('Anomaly after CO2 doubling')
for ax in [ax1, ax2]:
ax.set_xlim(-90, 90)
ax.grid()
ax.legend()
Seasonally, the dominant balance by far is between solar radiation and heat storage!
These notes largely follow Chapter 4 of Hartmann (1994) "Global Physical Climatology", Academic Press.
Turbulent fluxes of heat: eddy fluxes of heat and moisture at some level in the atmospheric boundary layer
$$ \text{SH} = c_p ~\rho ~ \overline{w^\prime T^\prime} $$$$ \text{LE} = L ~\rho ~\overline{w^\prime q^\prime} $$where $c_p$ is the specific heat of air at constant pressure, $L$ is the latent heat of vaporization, $\text{SH}$ is the sensible heat flux and $\text{LE}$ is the latent heat flux.
From theory of boundary layer turbulence, we suppose that the eddy heat fluxes is related to boundary layer temperature gradients, as well as the mean wind speed:
$$ \text{SH} = c_p ~\rho ~ C_D ~ U \left( T_s - T_a \right) $$where $T_s$ is the surface temperature and $T_a$ is the air temperature at some reference height above the surface. $U$ is the wind speed at the reference height, and $C_D$ is a dimensionless aerodynamic drag coefficient.
$C_D$ will depend, among other things, on the roughness of the surface.
Similarly, we assume that the latent heat flux is related to boundary layer moisture gradients:
$$ \text{LE} = L ~\rho ~ C_D ~ U \left( q_s - q_a \right) $$where $q_s$ is the specific humidity of air immediately above the surface, and $q_a$ is the specific humidity at the reference height.
In general the transfer coefficients $C_D$ could be different for sensible and latent heat flux, but empirically they are found to be very similar to each other. We will assume they are equal here.
The Bowen ratio is a dimensionless number defined as
$$ B_o = \frac{\text{SH}}{\text{LE}} $$i.e. the ratio of sensible heat loss to evaporative cooling.
From the above plots, the Bowen ratio tends to be small in the low latitudes.
Over a water surface or a very wet land surface, we may assume that the mixing ratio of water vapor at the surface is equal to the saturation mixing ratio $q^*$ at the temperature of the surface:
$$ q_s = q^*(T_s) $$Recall that the saturation vapor pressure $q^*$ is a sensitive function of temperature through the Clausius-Claperyon relation. (It also depends on pressure)
Let's approximate the mixing ratio for saturated air at the reference height through a first-order Taylor series expansion:
$$ q_a^* \approx q_s^*(T_s) + \frac{\partial q^*}{\partial T} \left( T_a - T_s \right) $$The actual mixing ratio at the reference height can be expressed as
$$ q_a = r ~ q_a^* $$where $r$ is the relative humidity at that level.
Then we have an appoximation for $q_a$ in terms of temperature gradients:
$$ q_a \approx r \left( q_s^*(T_s) + \frac{\partial q^*}{\partial T} \left( T_a - T_s \right) \right) $$Substituting this into the bulk formula for latent heat flux, we get
$$ \text{LE} \approx L ~\rho ~ C_D ~ U \left( q_s^* - r \left( q_s^* + \frac{\partial q^*}{\partial T} \left( T_a - T_s \right) \right) \right) $$or, rearranging a bit,
$$ \text{LE} \approx L ~\rho ~ C_D ~ U \left( (1-r) ~ q_s^* + r \frac{\partial q^*}{\partial T} \left( T_s - T_a \right) \right) $$The Bowen ratio is thus
$$ B_o = \frac{c_p}{ L \left( \frac{(1-r)}{\left( T_s - T_a \right)} q_s^* + r \frac{\partial q^*}{\partial T} \right)} $$Notice that if the boundary layer air is saturated, then $r=1$ and the Bowen ratio takes on a special value
$$ B_e = \frac{c_p}{ L \frac{\partial q^*}{\partial T} } $$When the surface and the air at the reference level are saturated, the Bowen ratio approaches the value $B_e$, which is called the equilibrium Bowen ratio. We presume that the flux of moisture from the boundary layer to the free atmosphere is sufficient to just balance the upward flux of moisture from the surface so that the humidity at the reference height is in equilibrium at the saturation value.
Recall that from the Clausius-Claperyon relation, the rate of change of the saturation mixing ratio is itself a strong function of temperature:
$$ \frac{\partial q^*}{\partial T} = q^*(T) \frac{L}{R_v ~ T^2} $$Here the quasi-exponential dependence of $q^*$ on $T$ far outweighs the inverse square dependence, so the equilibrium Bowen ratio decreases roughly exponentially with temperature.
The following code reproduces Figure 4.10 of Hartmann (1994).
from climlab.utils.thermo import qsat
T = np.linspace(-40, 40) + const.tempCtoK
qstar = qsat(T, const.ps) # in kg / kg
def Be(T):
qstar = qsat(T, const.ps) # in kg / kg
dqstardT = qstar * const.Lhvap / const.Rv / T**2
return const.cp / const.Lhvap / dqstardT
fig = plt.figure()
ax = fig.add_subplot(111)
ax.semilogy(T + const.tempKtoC, qstar*1000, label='$q^*$')
ax.semilogy(T + const.tempKtoC, Be(T), label='$B_e$')
ax.grid()
ax.set_xlabel('Temperature (degC)')
ax.legend(loc='upper center')
ax.set_title('Saturation specific humidity (g/kg) and equilibrium Bowen ratio')
Bo_control = surface_budget['control']['SHF']/surface_budget['control']['LHF']
Be_control = np.mean(Be(runs['control'].variables['TS'][:]), axis=0)
fig = plt.figure(figsize=(16,10))
ax1 = fig.add_subplot(2, 2, 1)
cax1 = ax1.pcolormesh(lon, lat, Bo_control,
vmin=0., vmax=5. )
fig.colorbar(cax1)
ax1.set_title('$B_o$ (CESM control)', fontsize=20)
ax2 = fig.add_subplot(2, 2, 2)
cax2 = ax2.pcolormesh(lon, lat, Be_control,
vmin=0., vmax=5. )
fig.colorbar(cax2)
ax2.set_title('$B_e$ (CESM control)', fontsize=20)
ax3 = fig.add_subplot(2, 2, 3)
cax3 = ax3.pcolormesh(lon, lat, Bo_control - Be_control,
cmap='seismic', vmin=-10., vmax=10. )
fig.colorbar(cax3)
ax3.set_title('$B_o - B_e$ (CESM control)', fontsize=20)
for ax in [ax1, ax2, ax3]:
ax.set_xlim(0, 360)
ax.set_ylim(-90, 90)
ax.contour( lon, lat, topo.variables['LANDFRAC'][:], [0.5,0.5], colors='k');
On the difference plot, the blue colors indicate the actual Bowen ratio is smaller than the equilibrium Bowen ratio. This will typically occur for wet surfaces with undersaturated air.
The red colors indicate the actual Bowen ratio is larger than the equilibrium Bowen ratio. This typically occurs for dry surfaces where there is not enough water available to satisfy the energetic demand for evaporation.
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences, offered in Spring 2015.
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%version_information numpy, climlab